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COMPUTATIONAL ALTERNATIVES TO OBTAIN 


TIME OPTIMAL JET ENGINE CONTROL 
Abstract 


This work presents two computational methods to 
determine an open loop time optimal control sequence for 
a simple single-spool turbojet engine described by a set of 
non- linear differential equations. Both methods are 
modifications of TNridely accepted algorithms which can solve 
fixed time unconstrained optimal control problems with a free 
right end. Constrained problems to be considered have fixed 
right ends and free time. 

t 

Dynamic Programming, originally formulated by Bellman, 
is defined on a standard problem and it yields a successive 
approximation solution to the time optimal problem of interest. 

A feedback control law is obtained and it is then used to 
determine the corresponding open loop control sequence. 

The Fletcher-Reeves Conjugate Gradient Method has been 
selected for adaptation to solve a non-linear optimal control 
problem with state variable and control constraints. It 
uses gradient information to improve the performance index 
of a nominal trajectory toward the minimum value. The control 
sequence which produces this minimum trajectory is the open 
loop time optimal control sequence. 

The two above methods are theoretically and computationally 
extended to include the free time, fixed right end time optimal 
control problem with state variable and control constraints. 





Computer software is developed which can solve a general 
class of constrained non-linear time optimal control 
problems. It is shown that the two methods produce, similar 
solutions to the turbojet problem.- 
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CHAPTER I 


original 
OF POOR QUALITY 


PROBLEM DEFINITION 


1.1 Introduction 

This chapter describes the unconstrained fixed time 

optimal control problem which well accepted computational 

methods are able to solve. It is seen that the jet engine 

presents a tougher problem due to state, and control variable 

constraints and due also to the fixed right end which is the 

design objective. Initially an optimal control problem is 

carefully defined and then extended so that the jet engine 

problem is included in the new class of problems for which 

the computational methods are to be adapted. Although the 

jet engine is described by a continuous time model, discrete 

* 

time systems are studied here because the methods are 
developed to find time optimal open loop control sequences 
on a digital computer. 

3 .2 The Standard Discrete Optimal Control Problem 
Consider the n order, time invariant discrete time 

system 


x(t+l) = x ( t ) + f (x(t) ,u(t) ) (1.2-1) 

with starting time k and terminal time N. In the system 
with m controls, x(t) is an n-dimensional state variable 
vector and u(t) is an m-dimensional control vector defined, 
at each sampling instant. In general, various starting times 
and states are considered, but we always denote 

x.(k) = x (1.2-2) 

as our initial time and state of interest. The terminal time 
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may be fixed or may be defined as the first instant at which 
the system state reaches a designated target set S. A 
performance index 

N-l 

J,(x-,u) = Y(x(N)) + £ L(x(t),u(t),t) (1.2-3) 

K t=k 

t = k,k+l, . . . ,N-1 

is to be minimized with u(t) fU, a control set, and u is 
the control sequence 

u = u(k) ,u(k+l) , . . . ,u(N-l) (1,2-4) 

It is understood that in the function L(x( t) ,u( t) , t) , 
x(t) is dependent upon the u(t) choice. 

■ Although there are proven methods which can solve the 
above problem, we actually -wish to solve a class of problems 
which include free time and fixed right' end conditions in 
addition to state and control ‘constraints of the form 

A x(t-) + B u(t) — C (1.2-5) 

It is this class in which the jet engine control problem lies. 
1.3 Jet Engine Control Problem 

From an accurate description of the Drone engine in Cl} > 
Brennan in £ 2 ] has derived a seventh order model and further 
reduced it to a third order model for simulation on a TR-48 
analogue computer. The discrete time version, obtained from 
the continuous time equations by Euler integration, is shown 
below and the physical representations of the variables are 
listed in figure 1.1. In this model, P^(t+<at) is understood 
to be P^(t+l) and P^(t) is understood to be P^. 
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PHYSICAL REPRESENTATION OF VARIABLES 
N = rotational speed 
P b = burner density 
P 4 = burner pressure 
w f = fuel flow 

= compressor discharge mass flow 
T^ = compressor discharge temperature 
Figure 1.1 


P E 

P 4 (t+1) = P 4 + At{(. 93586 |r + 31^86)w f 

P 2 

- 53-86 ) 

b 

< 

P b (t+l) = P b + A t(3?.78w 3 - 38.^8P 4 + . 

2 

N'( t+1 ) = N + a t(1.258/N)( - ^J 2 ) 


+ 21 .- 4 - 35 ^ 1 ^ 


(1.3-1) 

6684-9& f ) 

( 1 . 3 - 2 ) 

(1.3-3) 


= I. 3009 N - , 139825 P 4 - 

.13982 J? 4 2 + .41688N 2 - .0899P 4 N ' (i .3— 4*) 

T 3 ~ .64-212 + .35788N 2 (1.3-5) 


rhe two constraints are 1, the surge margin constraint 


P 4 * 1.25N 


(1.3-6) 


and 2, the turbine inlet temperature constraint , 

P 4 £ 1.25P b (1.3-7) 

The variables in the above equations are normalized in 
LI about an equilibrium point (P 4 -=1,,-P b g =l t N g =l ) such that 
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f(P, ,P, , N ) = 0. The control objective is to find a 

discrete open loop fuel sequence which takes the system from 
an initial state x of windmill »774,N=.546l) 

to a final state x p of military thrust (P^=P,P- b =l ,N=l) in 
minimum time such that the surge margain and turbine inlet 
temperature constraints are satisfied. 

1.4 Reduction to Second Order Model 

The effect of w f on the system occurs primarily in 
equation (1.3-1); its effect on equation (1.3-2) is minor. 

In addition the main single influence on equation (1.3-2) 
is so the assumption that P^ can be controlled almost 
directly by ^ is made. Therefore, a reduction of the system 
to a second order problem is made and it is. expressed, in 
the discrete state variable form of section 1.2, as 

x 1 (t+l) = x 1 + &t-(.3?.78w 3 - 38.448u 1 + .66849) 

(1.4-1) 

2 ? 

X 2(t+l) = x 2 + *t(1.258/x 2 )( ^1_ - w 3 x 2 ) (1.4-2) 

X 1 

ft 3 = 1.3009x2 - .139825u 1 - 

— ■■ 

u 1 2 + .4l688x 2 2 - .0899^^2 (1.4-3) 

with the control constraints 


U^ £ 1.25 X 2 

(1.4-4) 

U 1 - 1 - 25 X 1 

(1.4-5) 


In this second order model, P^ is now understood to be the 
control variable with P^ = x^ and N = x 2 the state variables. 
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Equation (!.*}-< T) includes the constant term because at the 
equilibrium condition where the constant term would have 
its largest effect, should have a value close to one. 

The state variable constraints’ remain identical to the third 
order model constraints but with the following subtle 
difference; now, the control is constrained by its position 
in the state -space. 

1 .5' Scope 

Three accepted methods which can solve the unconstrained, 
fixed time,- free right end problem described in section 1.2 
are Dynamic Programming, the Discrete Minimum Principle, 
and ordinary mathematical programming which minimizes a 
function of many variables. In this thesis, two of these 
methods are adapted to solve the general class of fixed right 
- end, time optimal control problems with state variable and 
control constraints. 

t 

In Chapter II, a successive approximations technique 
extends Dynamic Programming to produce a time optimal feedback 
control law from which an open loop control sequence can be 
determined. In Chapter III, non-linear programming is used 
to solve the jet engine problem and produce an open loop control 
sequence. Constraints are difficult to handle in the 
Conjugate-Gradient approach, but a simple feasible directions 
technique is used in conjunction with penalty functions to 
give satisfactory convergence without jamming. Careful 
inspection of the respective solutions shows that each method 
yields the same time optimal control sequence. 
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CHAPTER II 
DYNAMIC PROGRAMMING 

■ 2.1 Motivation 

Dynamic Programming is a method well suited- for solving 
low order, fixed time, free right end optimal control problems 
with state and control constraints. The feedback control 
solution and the cost function are calculated in successive 
steps as the state and control constraints actually reduce 
the required number of computations per step. The main 
difficulty is to include the fixed right end, free time 
problems into the class of problems which Dynamic Programming 
can solve. The appropriate modification of the theory is 
given below- 

Pontryagin in £ 3*} has shown that time optimal control 
problems very often tend to have bang-bang solutions. The 
Dynamic Programming method is also well suited to the bang- 
bang solution because the control is quantized and each 
control candidate is tested at each point in the state space 
to determine the optimal control. Therefore, regions of 
similar controls are expected and boundaries separating these 
regions are not surprising. Suggestions on accurately 
representing these boundaries are presented later in this 
chapter. 

2.2 Theory 

Let K denote the integers, R the real numbers, X an 
aroitrary nonempty state set, and S^rXxK an arbitrary nonempty 
target set. Suppose we also have a control constraint U(x) set 



7 


which depends on the state x and is such that controls in 
U(x) guarantee that the next state is in the state set X. 

Consider the discrete time system of section 1.2 and 
the performance index J^(x,u) in (1.2-3’) to be minimized. 

We assume that 

V k (x) = J k (x,u) (2.2-1) 

exists, where u denotes any admissible sequence u(k),u(k+l), 

. . . ,u(N-l) with values u(t) 6 U(x( t) ) . Only the above 
assumption is required to prove Bellman's Conditions [4l , £ 5^ , 
and [6j 

V k^ = u^UCx)^ L ^ x,u,k ^ + V k+l^ x+f ^ X,U ^ 1 

(x,k)^S (2.2-2) 


V k (x) = X(x) 


(x,k) <?■ S 

which are necessary and sufficient for optimality. These 
are the usual equations of Dynamic Programming. A particular 
u = v(x,k) <J U(x) which minimizes the expression above can 
be used to make up v(x,k) for all (x,k) £ S and comprises an 
optimal feedback control law. 

The optimal open .loop control sequence u can be determined 
from v(x,k) as follows. Let x y (t) be the optimal trajectory 
obtained using v(x,k) as a feedback control law. Then 

u(t) =v(x v (t),t) (2.2-3) 


and 


u = u(k) ,u('k+l ) , . . . ,u(N-l) 


(2.2-4) 


An alternative method for solving problems of this class 



8 


called, "successive approximations", was suggested by Bellman 
and later developed by Leake, Liu, and Richardson in [6] 
and \ 7 ~\ • 

Let V^ n (.x) be any function such that V^ n '(x) - V k (x) an< ^ 
let v n (x,k) be a control law which results when performing 
the minimization 

uVu<x) [ L(x » u ,k) + V k+1 n (x+f(x,u))] (2.2-5) 

(x,k ) 4 S 

Then it is shown in that if V k n+ ^(x) is the performance 
function resulting from v n (x,k), we have 

V k (x) * V k n+1 (x) ± \ n (x) (2.2-6) 

and further V k (x) converges monotonically to V k (x) in a 
finite number of steps although -each (x,k) may require a 
different number of steps. 

In the special case where V k (x) is independent of k, 

V k (x) = V(x), the approximating sequence V k n (x) may be taken 
as independent of k also to yield V k n (x) = V n (x) and (2.2-5) 
becomes 

V n+1 (x) = ^ i J 1 U ( x ) j^L(x,u,k) + V n (x+f (x,u) ) 

(x,k)/ S (2.2-7) 

Problem of Interest 

Recall now that in the jet engine problem (sections 1.3 
and 1.4) we wish to find the minimum number of steps (and 
associated optimal controls) that it takes to drive the system 
x(t+l) = x ( t ) + f(x(t) ,u(t) ) 


( 2 . 2 - 8 ) 
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from x(k) = x to a fixed state x p with associated state and 
control constraints x(t)£ X ana u(t) £ U(x( t) ) . Note that the 
minimum number of steps to reach the target state x p does 
not depend on the starting time k. Thus the minimum number 
of steps can be expressed as ‘ • • 

V k (x) = V (x) (2.2-9) 

a function depending only on x. Similarly the control law 
v(x,k) = v(x) . 

Ordinary Dynamic Dynamic Programming Problem 

We now define a standard Dynamic Programming problem 
which yields a successive approximation solution to the time 
optimal problem of interest. Simply let the system be the 
same as before with target set 

S = S 1 Vj S 2 (2.2-10) 

where S is illustrated in figure 2.1 and 

51 = { (x,k) : k = 0, xtX j (2.2-11) 

5 2 ““ , f(x,k) : k ^ 0, x = Xp^ 

Let X and U(x) be as above, let N=0, let L(x,u,k) = 1 to 
provide a time penalty to the system, andlet X(x) be chosen 
such that 

Y(x) * V (x) (2.2-12) 

X(x p ) = 0 

Then Bellman's conditions are as in equation (2.2-2) and 
the problem can be solved as an ordinary Dynamic Programming 
problem. Keep in mind however, that (x p ,k)£S so we always 
have 

ORIGINAL PAGE IS^. 

OF POOR QUALITY" 
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TARGET SET FOR SUCCESSIVE APPROXIMATIONS DEVELOPMENT 



-3 -2 -1 0 


Vertical lines represent the set X 
defined at each k. 


Figure 2.1 
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V k (x p ) = 0 (2.2-13) 

for any k - N = 0. 

Successive Approximation Solution 

To show that the ordinary Dynamic Programming solution 
above amounts to a successive approximation solution to the 
problem of interest, note that the actual target set of 
interest is simply the set of all (x,k) where x = x p . We 
can establish a successive approximation of V(x) by choosing 

V°(x) = Y(x) (2.2-14) 

Then, if we equate 

V k (x) = v~ k (x) , k = 0,-1, -2,-3, .. . (2.2-15) 

we see that Bellman's Conditions are equivalent to 

V n+ 1. = [" L (x,u,k) + V n (x+f (x,u) ) j (2.2-16) 

(x,k) ( S 

n = 0,1,2, o . . 

which is simply the successive approximation equation. This 
establishes that the Dynamic Programming problem actually 
gives a successive approximation solution, with 

V n (x)J^V(x) (2.2-17) 

and associated feedback control law 

v n (x) — ^v(x) (2.2-18) 

In practice, the state set is discretized and inter- 
polation is used to get approximate solutions, but convergence 
still occurs. Choosing some k T as a time when convergence is 
adequate, we equate 
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V(x) 


(2.2-19) 


where V(x) represents the number of time steps needed for 

v(x,k r p)=v(x) (2.2-20) 


to take the system from initial state x to final state Xp. 
V(x) also reveals those initial states x for which a control 
sequence such that u(t)£ U(x(t)) cannot be found to move the 
system to ‘the final s-tate Xp. The uncontrollable regions 
are simply found by observing that U(x) = 0 f or the states 
in those regions. In the computer program, a penalty cost 
is arbitrarily assigned to those states and a suitable 
representation is given to the corresponding feedback control 
laws to indicate that no admissible control exists. 

2.3 State Variable Quantization 

The first computational requirement in using a dynamic 
program is to determine an adequate quantization, independent 
of k, of the state variables and the control variables. In 
the jet engine problem, the time constants of P fe and N, 
obtained from linearizing the third order model about the 
design objective, (see [23 or section 3«5) differ by an 
order of magnitude. This difference implies that will be 
able to react more rapidly than N, To complicate matters, a 
suitable choice for the time step size, for the Euler 
integration must be -made to restrict the motion of one time 
step to one increment in each state variable. 

Larson in C5] presents a relation which is useful for 
determining a proper state variable quantization. If we let 
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,F(x( t) ,u( t ) ) = f(x(t),u(t))/^t 


(2.3-1) 


and let x. (t) represent the i^ component of x(t), let u.(t) 
represent the j th component of u(t), and let represent 
the i th discrete time function, then Larson's relation is 


. mm 

at = , 


x.(t) 




(2.3-2) 


The above equation requires that F^(x. (t) ,u^ ( t) ) be scanned 
over all possible states and controls to determine the maximum 
motion of each state variable per time step. A short 


= P. 


b max 


= 40 


computer program easily determines that F^ max 
and F 0 = N = 2 for the jet engine problem. From this 
information, the values of ^t, ^P-^, an< ^ are ( ^ e '* :erm lned . 

In problems with state variable and control constraints, 
each possible control must be in the set U(x) before it is 
tested to find the best u(t). The quantization of these 


problems should include a large number of state variable points 
which land on the boundary of the constraints- and it is 
necessary that those points in be members of the quantized 
grid. In the jet engine problem, the discrete quantization 
must include the point (P- b =l,N-l) and several "points which 
satisfy constraints (1.4-4) and (1.4-5) with equality. 

2.4 Initial Cost Function 

The initial assigned cost function, V^(x), is designed 
to force the system to rapidly approach Xp. One approach is 
to assign an arbitrary large constant to all x / x ? but it 
is found to decrease the convergence rate of the successive 
approximations solution because the large initial constant 



introduces significant interpolation errors into the algorithm 
If V(x) denotes the minimum number of time steps required 
for the feedback control law to take the system from an 
initial state x to a final state x F , equation (2.2-12) 
establishes a lower bound for the initial constant which is 
the maximum time expected divided by the time step size. A 
lower constant function increases the convergence rate of 
the successive approximations method and also reduces the 
errors induced by interpolation. The lower bound can be 
found with a few trial Dynamic Programs or with preliminary 
equation study. For the jet engine problem, the optimal time 
required to move from windmill to military thrust is slightly 
less than one second, or approximately .8 seconds. 

Consequently the constant is set equal to 1/^t. 

2.5 Interpolation 

Recall that V (x) will have numerical representations 
stored on the computer only for those states x which lie on 
the discrete quantization grid. If a particular control 
u£U(x) forces the motion 

x' = x + f(x,u) (2.5-1) 

to a state x 1 not defined on the discrete grid, interpolation 
is required to approximate V (x 1 ). Figure 2.2 illustrates 
the general case and presents an interpolation formula which 
produces good results for the jet engine problem. 

Early results have indicated a tendency for the cost 
function associated with the members of V - ^(xj,) = V^(x F ,k) 
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i, 3+1 


v i+l,j+l 

b 

V 

P 



t 


V. • 
J 

' 1 1 ■ 1 

a 

V i+1, j 


Vp = (l-a)(l-b)V i>3 + b(l-a)V lj3+1 + a(l-b)V 1+1> . + abV i+lj . +1 

Figure 2.2 


to incur with each iteration small positive values which 
propagate through the algorithm and grow larger. These 
errors arise from quantizing and interpolating non- linear 
system equations and cost functions. To minimize this effect, 
V n (xp) is reset to zero after each iteration. 

In figure 2»3> a flow chart, describing the general 
Dynamic Programming Method with the successive approximations 
extension, is presented. The actual computer program, written 
in FORTRAN IV for use on the FORTGI compiler in Notre Dame’s 
IBM 370/158 computer system, is listed and documented in 
appendix A. 

2.6 Control Regions for the Turbojet Engine 

Recall the reduced jet engine problem of section 1,^ 



Plow Chart for Dynamic Programming Method 


Initialize VOID 
Enter state grid parameters,, 
time step size, and number 
of iterations 


Iteration Index 
State Space Index 



A 

/ 


NOT < 0 
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Stop' 


Figure 2.3-A 
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Figure 2«3-B 
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which has one control variable P^ (or u^ ) and two state 
variables, P^ (or x^) and N (or x^). The control set U 
consists of real numbers which lie in the interval (0„5,' 

1.25) and various state sets X are considered which have 
members surrounding the design objective. Recall that the 
quantization is independent of k so that the design objective 
can loosely be referred to as the target. As the Dynamic 
Programming Method is run, the successive approximation 
solutions show that the controllable region of P^ grows 
quickly while the controllable region of N grows slowly. 

This is due to the vast difference in the time constants of 
P^ (Tp =. 0309 ) and N (T^=.352). A state set X consisting of 
a narrower range of N (0. 8-1.1) and a wide range of P fe (0.5- 
1.8) produces the most meaningful results. For greater 
detail near control region boundaries, smaller quantization 
increments and state sets X defining smaller regions can be 
studied. They yield more information about the boundaries 
and the feedback control regions are more clearly defined. 

The time optimal feedback control solution is illustrated 
on a state variable map in which three distinct control regions 
can be seen in figure 2.4. The sharp boundaries which are 
examined closely in section 2.7 and which separate these 
regions are indicated by solid lines while the motion of the 
system, x v (x,k,p), due to the feedback control law, v(x,k,p), 
is indicated' by dotted lines for several starting states x. 

The design objective Xp is denoted by a circle and every tenth 
of a second along the trajectories is marked by crosses. 



State Variable Motion when Time Optimal Control is Appli 



Figure 2«4-A 
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In region 1, the control variable rides the surge margain 
constraint. This causes a rapid decrease in P^ with a little 
increase in N. When the motion reaches the 1:2 boundary, the 
control switches to ride the turbine inlet temperature 
constraint'. In region 2, all starting states have motions 
which take the system to a common main path which does not 
lie on a region boundary. Along this path, N increases to 
its design objective but only a small increase in is 
observed. -At a point determined by the 2:3 boundary, the 
control jumps to its minimum value and the system motion 
rapidly approaches the target where the control assumes a 
value' of one. It is- this rapid switching action which 
Pontryagin refers to a's a bang-bang control. 

i 

Region y is -the area in which the value of N is higher 
than the military thrust value. Physically, this area is of 
little interest, but it reveals an interesting control pattern 
which has no obvious analytic properties such as those in 
regions 1 and 2. The results in figure 2.4-B are obtained 
from a control set U 6 ( *75-1*25) » and state variable set 
N £ (.97-1*025) and P^^- (.6-1.5)* The time optimal feedback 
solution reduces N rapidly by reducing the control variable 
to its minimum value., however the low control simultaneously 
increases the value of P^ above the design, objective. The 
dotted line contours indicate when the control value is 
continuously raised to cause the motion of the system to 
move along the 1:,3 boundary, where the control follov^s- the 
surge margin, directly to the target. 
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2 . 7 Accurate Boundary Descriptions 

At certain sections in the state map, the time optimal 
feedback control solution reveals the location of a boundary 
where the trend of controls changes. Between regions 1 and 
3, a line is implied for which states above have an optimal 
feedback control value v(x) =: 1 . 25 N and for which states 
below have feedback control values which satisfy the surge 
margain constraint with equality. The exact boundary is the 
trajectory which, formed by the motion of the system due to 
controls lying on the surge margain constraint, proceeds 
directly to the target. This is determined in figure 2.5 by 
integrating the system forward with.^a few selected states as 
initial conditions. In the neighborhood of this boundary, 
the control contour is continuous, or in other words, there 
is a smooth transition of feedback control values from region 
3 into the boundary. 

The 2:3 boundary separates two regions having radically 
different optimal control strategies. States x which lie 
in region 2 have v(x) = 1.25P b and states x which lie on and' : 
above this boundary in region 3 have corresponding minimum 
feedback control values. Figure 2.6 illustrates the resulting 
boundary as a function of m ^ n which can be chosen arbitrarily 
in the reduced jet engine problem. Its most realistic value 
depends upon how quickly P^ decreases when w f in the third 
order model is suddenly reduced to zero. 

The boundary separating regions 1 and 2 is simply the 


line P b « N. 
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2.8 Summary of Second Order Model Controls 

Figure 2.7 shows the optimal time plots of the reduced 
jet engine problem to which the time optimal feedback control 
laws obtained from the dynamic programming algorithm are 
applied. P^ has three control rules. Initially, P^ follows 
the surge margin constraint until P-^ is lowered sufficiently 
to equal N, at which point P^ will follow the turbine inlet 
temperature constraint. Physically, these two rules represent 
the maximum allowable throttle without stalling or overheating 
the jet engine. At .79 seconds, N has a value slightly 
above the design speed but P fc will have only 80 % of its 
design value. The point is determined by the 2:3-boundary 
where P^ is dropped .t-o P^ min for an instant and then returned 
to a value of one. This' discontinuous rapid switching action 
is an example of a bang -bang control -and it -moves P^ directly 
to the design value in only .01 seconds. Therefore, the jet 
engine is controlled from windmill to military thrust in 
.8 seconds. 

2.9 Application to Third Order System 

From the reduced system information, the optimal time 
responses for the third order description of the Drone engine 
can be calculated. To determine the proper w^, it is required 
to find the desired P^ from the time optimal feedback control 
law and calculate the value of w^ which will force P^ to 
assume values which follow the time optimal control law. 

When the trajectory reaches the 2:3 boundary, w^. is suddenly 
reduced to a zero value which causes P^ to decrease and P^ to 
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increase to its design value. Then, must be a positive 
impulse to raise P^ to the target rapidly without changing 
N or P fe . Once this is accomplished, is set equal to one 
and the system is in the military thrust state® 

The time responses for the third order model are plotted 
in figure 2.8. The initial impulse of w f forces P^ up to 
the N constraint in one time step. 







29 


CHAPTER III 

THE MODIFIED FLETCHER-REEVES CONJUGATE GRADIENT METHOD 
3*1 Motivation 

Ordinary non-linear mathematical programming is a well 
known method which can solve an unconstrained optimal control 
problem with fixed time and a free right end by minimizing 
a performance index, J(u), which is a function of several 
variables. This method determines the optimal open loop 
control, sequence and its computational requirements increase 
only arithmetically wi.th the order of the system under study. 

In this chapter, a feasible directions idea is presented in 
conjunction with penalty functions to extend the general class 
of problems which the Fletcher-Reeves Conjugate Gradient Method 
can solve to that which includes a constrained time optimal 
. control problem with a fixed right end. The developed 
computer software is listed and documented in appendix B for 
use to determine open loop time optimal control sequences for 
the general class of non-linear free time, fixed right end 
optimal control problems with state variable and control 
constraints. 

3-2 The Fletcher - Reeves Conjugate Gradient Method 

Consider the n^* 1 order, time invariant discrete system 
having m controls 

x(t+l) = x{ t ) + f (x( t) ,.u(t) ) (3.2-1) 

and 

x(-0) = x (3.2-2) 

with starting time zero, terminal time N, and f(x(t),u(t)) 
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is continuously differentiable over the entire state space 
and control domain. The performance index to be minimized 
is 

N-l 

J { u ) = K(x(N)) + L(x(t) ,u(t) ,t) . (3*2-3) 

t=0 

t = 0,1,2, . . . ,N-1 

which is a function only of u given a starting state x where 
u denotes the control sequence 

u = u(0) ,u{l) , . . . ,u(N-l ) (3*2-^-) 

and 

u( j) = u x ( j) ,u 2 ( j) , . . . ,u m ( j) (3.2-5) 

so that u is a sequence of Nm real numbers. 

This general clas's of unc onstrainecl free right end fixed 
time optimal control problems can be solved by the Fletcher- 
Reeves method which uses a combination of steepest descent 
and conjugate gradient techniques presented in Us] to minimize 
a directionally convex multivariable function, J(u). In 
figure 3*1 » a flow chart for this algorithm is presented and 
the appropriate equations are described below. 

Suppose a nominal control sequence u with a resulting 
state variable trajectory and .performance index is selected. 

The gradient of the trajectory is calculated while the adjoint 
system equations are solved in reverse time according to 

y(t) = y(t+l) + y(t+l) V x f (x(t) ,u(t) ) + tf x L(x(t ) ,u(t) ) 


y(N) = V x K(x(N)) 


t = N-l , N-2 , . . . ,2,1 


( 3 * 2 - 6 ) 












32 


with y(t) an array of n length, f (x( t) ,u(t) ) is the Jacobian, 

A 


V f (x( t) ,u(t) ) 

A 


and 

£? x L(x(t) ,u(t) ) 
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(3*2-7) 


( 3 * 2 - 8 ) 


The gradient is the Nm length sequence 
£ 


jv (a) = t 7 u(0) J( li) 7 u(1) J(u) 


( 3 . 2 - 9 ) 


and each component can be calculated from the relation 

v u(t) J( -) = y(* + i ) (x(t) ,u(t) ) + v u L (x(t) »u(t) ) 

(3.2-10) 

with V u f (x( t) ,u(t) ) and V^L(x( t) ,u( t) ) similarly defined. 

It is incidental that the Fletcher-Reeves method under 
appropriate convexity assumptions actually solves, the discrete 
minimum principle for the above class of problems. If the 
Hamiltonian of the system is defined by 

H(x(t) ,u(t),y(t+l) ) = y( t+l)f (x(t) ,u(t) ) + L(x(t),u(t)) 

( 3 . 2 - 11 ) 

then 

V u ( t )J(u) = \7 u H(x(t) ,u(t) ,y(t+l)) (3.2-12) 

which is zero along the optimal trajectory, corresponding to 



the Hamiltonian "being minimized. 

During each first inner iteration, a line search is 
performed along the conjugate direction of the gradient 
cL - -£ to find the best d such that 

h(* *). = “i n 0 J(u + <* d) (3.2-13) 

For each successive inner iteration, the direction of the 
line search is 

d* = -g* + Ed (3.2-14) 

a linear, combination of the present gradient and previous 
direction with 

£ . 

B = (3.2-15) 

■K E 

and g* is the present ' gradient , g the previous gradient, and 
d is the previous direction. -After N inner - iterations, the 
direction is again set equal to the conjugate direction of 
the gradient and the present performance index is compared to 
the previous one to determine if another set of inner 
iterations is required. If, during a set of inner iterations, 
the line search produces ,an.^ * .= 0-, the .performance index 
can not be improved with further inner iterations so the 
algorithm skips to the next outer iteration. The control 
sequence u which minimizes J(u) is the optimal open loop 
control sequence. 

3.3 The Extended General Class of Problems 

We actually wish to solve a constrained time optimal 
control problem with -a fixed right -end . -Recall the reduced 



jet engine problem of section 1.4 with control constraints 
and a design objective x F = (P tj =l,N=l). A combination of 
a simple feasible directions idea used in conjunction with 
penalty functions will extend the capabilities of the Fletcher 
Beeves' Conjugate Gradient Method to solve the general 
constrained time optimal control problem with a fixed right 
end. 

3*4 Feasible Directions Modification 

Several barrier functions were tried which would heavily 
penalize the system if a control constraint is violated. All 
attempts to find a continuously differentiable function to 
force the control variable, P^ t to obey its constraints failed 
■so great is' -the tendency of -the '-system to overshoot them. The 
overshoot path always yields a lower performance index than 
any legal control .sequence would yield. As the barrier 
functions are steepened, computer overflows appear when the 
performance index is differentiated, and therefore barrier 
functions do not provide a solution to the overshoot problem 
in the reduced jet engine model. 

Hence, barrier functions .are abandoned and a simple 
feasible directions idea is added to the algorithm. If a 
member of a control sequence violates its constraints, the 
control is reset at the nearest constraint boundary, 
making it impossible for any control member to violate the 
control constraints. This subroutine is called from all 
trajectory calculations, and especially from those in the 
line search. -The- -'restricted ^algorithm -still uses gradient 
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Information to compute the steepest direction, but it restricts 
the magnitude of that direction which can be added to each 
member of the nominal control sequence. It should be noted 
that the gradients are no longer zero along the optimal 
path and also that a jamming possibility -exists. 

Recall in the reduced jet engine problem' that there 
are two oontrol constraints. In addition we constrain the 
value of P^ to be at least .5 to prevent negative values 
which would increase P^ and N arbitrarily. This constraint 
provides an equivalent problem to that solved by Dynamic 
Programming in chapter II. 

3*5 Performance Index Function 

A performance index function is used to convert the 
Fletcher-Reeves Method to solve a free time, fixed right end 
optimal control, problem. Since, in- the optimal time problem, 
the state variables are desired to reach x^, in minimum time, 
a performance index function, which penalizes the system for 
not being at x p , takes the form 

N-l n N-l , 

2 L(x(t),u(t)) = % 1 cAxAt) -,T.) 2 (3.5-1) 

t=0 j=l t=0 J J J 

where x F = (T^ ,T 2 , . . . ,T n ) and c = (c 1 ,c 2 , . . . ,c n ) are selected 

weighting constants. The squared criterion is chosen to 

provide a symmetric convex penalty function about the state 

x F which is continuously differentiable. Initially the c^'s 

were assigned the same value but this choice is found to 

produce an oscillating time optimal solution, a condition 
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not suggested .by the Dynamic Programming solution of chapter II. 

To remove the oscillation effect and obtain an accurate 
time optimal open loop control sequence, the c^s are 
chosen so they are in proportion, to the time constants of 
the state variables, which can be obtained by linearizing the 


system about x p . If we represent the A,B linear state 
description of small deviations about the design objective as 

-Mr * a <£ x + B £u (3*5-2) 

a x 

with 

A - (x p ,u F ) (3*5-3) 

then the time constants are the inverse of the real parts of 


the eigenvalues - of the Jacobian evaluated at x p and u p . In the 
jet engine problem, T N = .352',' T p = .0309, and T p = .0123 
and the resulting- c^ ratios become 1000 : 90 s 35 * The slowest 
state variable, N, should have the higher penalty because it 


requires the most time to react. In the reduced problem, 
is the control and is not assigned a penalty to find a 
time optimal control sequence for the problem. 

A prior knowledge of the approximate optimal time is 
needed to choose a good time step size and state variable 
array length. Chapter II has already shown that the jet 
engine can move from windmill to military thrust in .8 seconds 


so a t =s .005 and an N = 200 are used in the program. 

For the optimal time problem, the early states in the 
trajectory are given a small penalty with respect to that 
given to later states. This gives the system freedom to 



reach the design objective. If we let t be the time step of 
interest and let N time steps exist, then the performance 
index functions become 
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L(x(t) ,u(t) ) = N (™) 9 ^ c , (x . { t ) - T .)' 2 (3*5-4) 

j=l 3 2 J 


and 


n 

K(x(N) ) = N 5 c.(x,(N) 

j=l J J 



(3.5-5) 


/ 

A point 90 % along the trajectory is given a penalty only 38 % 
of that received by the endpoint. Computational experience 
has revealed that for problems in which the time optimal 
control sequence lies along the constraint boundaries for 
most of the state variable trajectory, a much lower exponent 
miJ.st.be used in (3*5-4) to provide higher earlier penalties. 

In run I of section 3 . 7 * an original exponent of 9 produces . . 
a solution in which the jet engine idles at windmill for most 
of the allowed time and only approaches the design value 
near the end. When the exponent is lowered to two, the 
solution is a control sequence which rapidly accelerates 
the engine. 

3*6 .Para me ter Sensitive Cubic Pit Line Search 
A quick, accurate line search to find 

= 2 l " 0 j(a +* 4) (3.6-1) 


is a vital component of the modified Fletcher-Reeves method. 
Lasdon in [ 9 ] has proposed a quadratic fit technique to find 
the best alpha. Also considered is a golden section line 
search but neither approach simultaneously allows a wide 
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range of possible d * values and requires relatively few 
trajectory calculations to determine a test point h(a). 

The line search proposed here fits a cubic polynomial 
through three test points and with derivative information 
algebraically solves for d * in one step* 

Consider the general situation pictured in figure 3*2 » 

Conditions 

1) h(a) ■£ h(0) 

2) h(b) 7 h(0) 

3) h E (0) * 0 

4) Oi a 4 b 

Figure 3.2 




To find A and B, note that 
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h(a) = h(o) + h*(0)a + Aa 2 + Ba 3 (3*6-6) 

h(b) = h(0) + h‘(0)b + Ab 2 + Bb 3 (3*6-7) 


and the expression for A is then obtained by multiplying 
(a/b) 3 to both sides of (3*6-7) and subtracting the result 
from (3*6-6). 


3 .3 

A — h(a) — h(b) + h(0) (—^— — 
b- 5 b J 

B - “ h(0) - h t (0)a - Aa 2 


3 

1) + h c (0H-^ - a) 

x b 

( 3 * 6 - 8 ) 

(3*6-9) 


a 3 - 

The flow chart for the parameter sensitive cubic fit 
line search is illustrated in figure 3*3* The components 
of the line search are ‘explained in the remainder of this 
section. Recall from the Fletcher-Reeves flow chart that 
h(,0) is already known and h*(0) can be expressed as g d. 
1) Selection of “a” Value 

As an initial guess for “a 11 , let 


a IU|| * 10 INI *=% a = iF*fd| (3*6-10) 

with l| || denoting the usual norm of a vector. If the 
resulting h(a) satisfies condition 1 of figure 3*2, let 
b = and proceed to find a suitable fi b fl value; if not, 
divide "a" by ten and try condition 1 again* Sometimes when 
J(u) is very close to the minimum value, H a" can get minutely 
small with no decrease in h(a)« If this happens, set ci * 
to zero and continue with the main program. 



Parameter Sensitive 
Cubic Fit Line Search 

















2) Selection of “b 1 * Value 


If the resulting h(b) satisfies condition 2 and has 
a reasonable value, use the cubic fit routine to find oC # 
with equations (3»6-5) , ( 3 . 6^*8) , and (3„6-9). If h(b) is 
less than both h(0) and h(a) , move “b" values into "a' 1 , let 
b = 2b and try condition 2 again. This step moves'- the value 
of "a” closer to the minimum point and improves the accuracy 
of the cubic fit. I-f condition 2 is not satisfied and h(b) 
is greater than h(a), let c* * = a. Experience indicates that 
in this case, the function h(c£ ) has the shape indicated in 
figure 3 .^. Since much computer time is spent to find a 



Figure 3 .^ 


"b" which will satisfy condition 2, and since both of the 
values of "a 11 and “b" are close to the minimum point, let 
= a. It must be noted that in the jet engine problem, 
this case only appears at the first iteration when the 
trajectory approaches the control constraints. After the 
cubic fit, there is a final test to check that h(e*.*) is 
actually lower than h{a), which will be the usual case. 

In the jet engine problem, this line search works well 
and always improves the ..performance index with -each iteration 
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There seem to be many tests in the line search of figure 3*3 » 
however the answer is yes to all tests in the usual case and 
only three or four trajectory calculations are required to 
find an accurate Another advantage is the wide range 

of values which * can assume* It ranges from 1.0 x to 
,0015 in the jet engine problem or eight orders of magnitud-e. 
3,7 Second Order Jet Engine Study 

In solving the jet engine problem and comparing the two 
computational methods of this thesis, trajectories resulting 
from two important sets of initial conditions are studied by 
the Modified Fletcher-Reeves Conjugate Gradient Method, From 
chapter II, it is already learned that the jet engine can be 
controlled from windmill to military' thrust in close to ,8 
seconds, so for this case, the state variable arrays will 
have a length of 200 and a time step size of .005* From 
figure 2.4-A, an initial state (P^=lo774,N-,9) is-seen to 
require about .3 seconds so a state variable array length of 
200 is used with a time step of ,002 for this problem. 

A poor nominal constant control of P^ = .6 which decreases 
the rotational speed and increases the burner density is 
used in the non-linear programs because this nominal control, 
unable to influence the solution, yields a true test of the 
feasible directions modification and free time, fixed right 
end extension. The number of inner iterations per outer 
iteration is equal to the array size of the variables. Outer 
iterations are run until the performance index no longer 
decreases. If the line search produces an c* * equal to zero. 
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no further improvement can be made until a new gradient 
direction is calculated* Therefore, significant c.p.u. time 
can be conserved by skipping the remaining inner iterations 
when anc£* equal to zero is produced. Table 3«1 tabulates 
the action of the performance index as outer iterations of 
the three significant programs are run. The vast improve- 
ment in the first iteration represents the rapid action away 
from the poor initial -control -sequence toward a much improved 
open loop control law-. Further iterations refine the optimal 
solution and only slightly lower the performance index. An 
illustration of this is the bang-bang portion of run III, in 
which the time plot of -the .control sequence assumes a 
rounded shape, in the -'earlier iterations and refines later 
to a rectangular : shape. • - -- - 

Sun I directly attacks the jet engine problem posed i-n - 
section 1.4 and produces very encouraging results. A close 
inspection of the time response in figure 3.5 reveals that 
the control sequence follows the control constraints through- 
out most of the trajectory. Near the end, the control is 
reduced and increased but not with the rapid switching action 
suggested by the dynamic programming solution. The author 
hypothesizes that the algorithm has jammed here because the 
optimal direction is calculated with no awareness of the 
control constraints and also the gradients along those 
constraints do not approach zero. However, a great part of 
the solution has been calculated with run I and it agrees 
with that in nhapt.er II. Figure 3*6 reveals that run 'II 



ITERATIONS REQUIRED AND RESULTING 
PERFORMANCE INDEX 


OUTER ITERATION # 


PERFORMANCE INDEX 


Run I x = (P b =l c 774 ,N=. 5461} t ,005 


0 

1 

2 

3 

4 

5 


2639895.00 
107942 „5'6 

107791.62 
107725.00 
107724.44 
no change 


Run II x = (P b =1.774,N=.9) t = .002 


0 

1 
2 

3 

4 


255667 3.00 

333.52 

328.07 

327.53 

no change 


Run III x = ( P b = . 7833,, N=. 9853 ) t = .001 


0 

1 

2 

3 

4 


88460.69 
101.49 
90.17 
89.22 
no change 


Table 3«1 
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Time Response of Reduced Jet Engine 


Obtained from Run I 



Figure 3.5 
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has produced a similar situation and has also defined a state 
variable path which joins that formed in run I. This oc- 
currence verifies the existence of the common trajectory of 
Dynamic Programming^ region 2 which is. se.en in figure 2c4~A. 

The initial conditions of run III are obtained by 
assuming that the portion of the solution, of which, the control 
values lie on the constraints, Is valid and by choosing a 
point near the end of the boundary as an initial condition for 
another run. The, new state variable array dimension is 100 
and the time step size is .001. The time response plot in 
figure 3*7 reveals a bang-bang control sequence which is 
faster than that in figure 3*5 and 3*6 and is similar to the 
solution of chapter II. It shows that P^ switches directly 
to the P^ constraint from the turbine inlet temperature 
constraint. This later sequence added to the control sequence 
obtained from the first part of run I produces an open loop 
control sequence which is most similar to that of chapter II. 

The methods and results of section 2.9 can again be 
employed to calculate the third order time optimal control 
sequence which is already graphed in figure 2.8. Since -the 
reduced solutions are identical, so are the third order solutions. 

In figure 3.8, the discontinuous control law for the 


two methods is shown 






Figure 3«8 
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Optimal Time Responses 
Dynamic Programming Method 



Figure 3 
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CHAPTER IV 
CONCLUSION 

Two well known computational methods, Dynamic Programming 
and the Fletcher-Beeves Conjugate Gradient Method, are 
successfully adapted in this thesis to solve the general 
class of fixed right end, free time optimal control problems 
with state variable and control constraints. Both methods 
have determined similar time optimal open loop control 
sequences for the discretized version of the reduced jet 
engine problem. These solutions have bang-bang features 
and also control values satisfying constraints with equality. 

A successive approximations technique extends the 
capabilities of Dynamic -Programming to -produce a time optimal 
feedback control law which has described, in the jet engine 
problem, several -regions in the. -state -space • map for which the 
controls follow a common rule. Separating these regions are 
distinct boundaries whose locations are approximately indicated 
on the discretized state variable grid. Unfortunately, the 
number of required computations increases geometrically with 
the order of the system and there is a problem with picturing 
a multi-dimensional feedback control law. 

A feasible directions idea, used in conjunction with 
penalty functions, extends the Fletcher-Reeves Conjugate 
Gradient Method to produce a time optimal open loop control 
sequence for a fixed right end optimal control problem with 
state variable and control constraints. The solution of the 
reduced jet -engine -problem again .shows the control values 
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lying on the control constraints for most of the state 
variable trajectory. Using a point near the end of the 
constraint region as- an initial condition, a discontinuous 
bang-bang control sequence, identical to that suggested by 
the feedback control law of Dynamic Programming, is obtained. 
This method is well suited to higher order problems because 
the number of computations required to solve higher order 
systems increases only arithmetically with- the system order. 

The computer software which produced the main results 
of this thesis has been generalized to solve any general 
problem in this class.. I-n the appendix, the software programs 
are listed and documented for general use to solve fixed 
right end, free time optimal control problems with state 
variable and control constraints,. 



APPENDIX A 


This appendix lists and documents the software which 
implements the successive approximations technique to 
Dynamic Programming. It is written with the reduced jet 
engine equations but it is carefully documented to describe 
a straightforward extension to the general constrained- 
time optimal control problem. The flow chart appears in 
figures '2.3 -A and 2.3-B. 
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DYNAMIC PROGRAMMING SOLUTION FOR SECOND ORDER SYSTEM 

THIS PROGRAM ALLOWS 'AM ARBITRARY SELECTION OR THE STATE SPACE 

TO BE STUDIEO 

DIMENSION VNEW ( 30 « 30.) »-V0LQ ( 30 , 30 ) ,U0PT(30,30) iALIST(30) 

NITT=DESlRED NO. OF ITERATIONS 
DELTA IS THE DESIRED TIME STEP- 


PB IS THE BURNER DENSITY STATE VARIABLE XI 
EN IS THE ROTOR SPEED STATE VARIABLE X2 

P4 IS. THE BURNER PRESSURE USED AS THE CONTROL IN THIS SYSTEM. 


FST IS THE FIRST POINT IN THE ALLOWABLE STATE SPACE 
LST IS THE LAST POINT IN THE ALLOWABLE STATE SPACE 
INC IS THE INCREMENT IN THE STATE SPACE • 

NOTE, RESTRAIN TO 30. -POINTS PER VARIABLE, 

the target points must be- on the grid 


ENFST=0.98 
ENLST=1 .05 
ENINC=0. 0025 
PBFST=0.b 
PBLST=l.b 
PBINC=0.05 
P4FST=0". 75- 
P4LST=1«25 
P4Il\iC = 0,125 
P4INC=P4'TNC/2. 

DELTA=.0.,.0.01 

NITT=T30 

ENTER THE .TARGET POINTS 

PBTAR=1. 

ENTAK=1. 

ER=l.«:E--6 

NEN'=1FIX,UENLST-ENFST+'EP>-/ENINC>+1 

NPB‘=IFIX (-(.RBLS-T-PBFST+EP ) /PB'INC ) +1 

NP4=IF1XL(P4L$T-P4FST+BP)/P4INC)+1 

NEN1=IFIX ( ( ENTAR-ENFST+EP ) /ENINC ) +1 

NPBl=IFIX( ( PBTAR-P8FST+EP ) /PBIWC ) +1 

NNEN=NEN+1 

NNPB=NPB+1 

NSP IS THE STARTING POINT FOR THE PB PRINTOUT 
NEP IS THE END POINT FOR THE PB PRINTOUT 


NSP=1 

NEP=NP3 

INITIALIZE THE PENALTY FUNCTION - HAKE AS LOW AS POSSIBLE 

DO 1'0 I=1,NNEN 
DO 10 J=1 , NNPB 
10 VOLDtl, J)~l ./DELTA 
VOLU ( NEN1 , - NPBl ) =0 • 

DO 15 1=1, NPB 

15 ALIST(I)=PBINC*(FLOAT(I-l) 1+PBFST+EP 


ITERATION INDEX 


DO 20 N=1,NITT 


STATE VARIABLE INDEX 


DO 30 1=1, NEN 

EN=ENINC* ( FLOAT ( 1-1) ) +ENFST 
DO 30’ J=1 , NPB 

PB=PBINC*(FLOAT< J-l) )+PBFST 
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NOT = 0 

VOTES=1./DELTA+10. 

CONTROL INDEX ■ 

DISCARD UNALLOWABLE CONTROLS 

DO AO K=l.iNP4 

P4=P4INC* < FLOAT ( K-l ).-) +P4FST 

THE CONTROL CONSTRAINTS ARE INDICATED HERE' 

IF (P4.GT •1.25001* EN) GO TO 40 
IF(P4,GT.1.25001*PB) GO TO- 40 

CALCULATE THE NEXT STATE ANO DISCARD UNALLOWABLE CONTROLS 

THE STATE EQUATIONS HERE ARE FOR THE REDUCED JET ENGINE PROBLEM 

Z=SQRT( R4*P4+0 .4T688*EN*EN-0 . 0 899*P4*EN ) 
W3DQT=1.30G9*EN~0.139825*P4-0. 13982*2 
PB,T=PB+DEL1 A*-('3'.7. 76*W3D0T-38. 44B*P4+0 . 66849 ) 

IF(PBT.GT.PBLST+EP) GO TO 40 
IF(PBT.LT.PBFST-EP) 60 TO 40 

ENT=EN+DEL1A*1.?SB/EN*(P4*P4/PB-W3D0T*EN*EN) 

IFTENT • GT , ENLST +EP ) GO TO 40 
IF( ENT. L r , ENFST-EP ) GO TO 40 
NOT=l\IOT+l 

INTERPOLATE TO FIND CORRESPONDING PENA7Y 

AI= ( FLOATT NE-N-1- ) *ENT-ENFST*FLOAT ( NEN ) +ENLST ) / ( ENLST-ENFST ) +EP 

bj= (Float ( npb-i j *pbt-pbfst*float.( npb i +pblst j / (PBlst-pbfst > +ep 

I COR- AT 

JC0R=BJ 

A=AI -FLOATT I.COR ) 

B=B J-'F.LQATT JCOR ) 

ONE=VOLD (JCOR, JC0R) 

TW0=VQLD (T-COK , JCOR+1 ) 

THREE,— VOLD (J-COR’.+fl,* JCOR:) 

F0UR=V0LDTIC0R+1, JCOR+1) 

V= < 1 . -A ) * < 1 . -B-) *0N£+ ( 1 . - A ) *B*TWO+A* ( 1 .--B ) *THREE'+ A*B*FOUR 
■DETERMINE IF THIS CONTROL IS BETTER' 

VTEST = 1-+V 

IFfVTEST.GT. VOTES) GO TO 40 
U0TES=P4 
VOTES=VTE'ST 
■40 CONTINUE 

FORM VNEW MATRIX FROM NEW DATA • 

IF ( NOT )'31 » 32 i 31 

31 IF.OVOTES. GT .l./DELTA-EP) GO TO 32 
UOPT ( I * J ) =U0.TES 

VNEW ( T» J) = V0TE'S 
GO . TO 30. 

NOTATION "FOR THE UNCONTROLLABLE STATE AND' FEEDBACK CONTROL LAW 

32 UOPT ( I i J ) =0 • 

VNEWUt J)=l. /DELTA 

30 CONTINUE 

RESET TARGET STATE PENALTY TO 0 
VNEW(NEN1,NPB1)=Q. 

MOVE VNEW INTO V0LD FOR NEXT ITERATION 

DO 50 1=1 i NEN 
DO 50 J=1 i NPB 
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50 VOLDU, J)=VNEW{Ii J) 
20 CONTINUE 

PRINT ROUTINE 


WRITE(6,100)NITT 

WRITE (6, 101) (ALIST(I) tI=NSPiNEP> 

WRITE(6<103) 

DO 33 K=1 i NEN 
KK=NEN+1-K 

XO=ENLST-ENINC*FLOAT ( K-l ) 
WRITE(6?l02)XDt (VOLD(KK,I ) , I=I\ISP , NEP ) 
WRITE (6. 101) (UOPT(KK.I) , I=NSP,NEP) 
WRITE < 6 * 103 ) 

33 CONTINUE 

WRITE (6, 101) ( ALIST(I) ,I=NSP.NEP) 

100 F0RMAT(5X. 'TIME STEP NUMBER* ..IT) 

101 FORMAT (6X*19F6»3) 

102 FORMAT ( IX i F5 « 3 1 19F6 . 1 ) 

103 FORMATtlX, ) 

STOP 

END 




APPENDIX B 


This appendix contains detailed software of the 
Modified Fletcher-Reeves Conjugate Gradient Method. A 
chart showing the relation of the main program with seven 
specialized subroutines is given below and the listings 
and documentations of the actual software are given on 
the following pages. 
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DIMENSION U(200,3) , US (200.3) »X< 201. 4) , XS(201,4) , 
?G(200,J3) , GS ( 200 , 3 )-i'D (200 *3) . ALIST-(>201 ) 

COMMON PCON ( 4 ) , -T ( 4 ) ', DELTA , KOUT * IND , NC , NX . NPA 


READ IN INITIAL DATA 

CARD 1 - NC = NUMBER OF CONTROLS 

NX = NUMBER OF STATE VARIABLES 
NPA = NUMBER OF POINTS PER ARRAY , 
NOUTsMAXIMUH NUMBER 'OF TOTAL ITERATIONS 
CARD 2 - INITIAL CONDITIONS FOR EACH STATE VARIABLE 
CARD 3 ■> PENALTY CONSTANT FOR EACH STATE VARIABLE 
CARD 4 ~ DESIGN OBJECTIVE FOR EACH STATE. VARIABLE 
CARD. 5 - TIME STEP SIZE 


READ C5 » # ) NC * NX , NPA , NOUT 
READ ( 5 , * ) (XU.-, J) tJ=l-,NX) 
READ (5,*) ( PCON ( J ) i J=1,NX) 
READ ( 5 , * ) ( T (U ) i J=1 , NX ) 
READ 05 1 * ) DELTA 
.DO 20 U=1 , l\IX 
20 XS(l-iU)=X(li J) 


INDEX IS THE .TOTAL ITERATION COUNTER 
KOUT IS THE' OUTER ITERATION' COUNTER 


• INDEX=0 
KOUT=0 
PTEST=1.E10 
EP=l.E-4 


ASSIGN A NOMINAL CONTROL 

DO 1 K=1,NC. 

DO 1 J=1-,NPA 

1 U ( J , K-) = 0 , .7 

2 CONTINUE 

IND IS 1 HE INNE-R- ITERATION COUNTER 


IND=0 

KOUT=KOUT+l 


•COMPUTE “THE STATE VARIABLE. TRAJECTORY AND THE PERFORMANCE INDEX 

CALL TRAJ( X , U ) 

CALL PERFM l X , U t PER") 

COMPUTE THE GRADIENTS 

CALL ADJNT<X,U.,G.) 

CALL PRIN(X,.U,G, PER', INDEX) 

DIRECTION IS THE NEGATIVE -.OF THE GRADIENT 


DO 3 K=1,NC 
DO 3 J=1 » NPA 

3 D< J , K ) =-G { J , K ) 

4 CONTINUE 
INDEX=INDEX+1 


LINE SEARCH TO FIND THE ALPHA .WHICH MINIMIZES 
H< ALPHA) = J(U:+ALPHA*D> 


CALL LNSCH(X«U,D,G, ALPHA, PER) 

IF ALPHA IS VERY SHALL, PROCEED TO THE NEXT CUTER ITERATION 

IF ( ALPHA. LT.l.E-12) INDEX=KOUT*NPA 
IF ( ALPHA, L'T , 1 . E-12 ) IND=NPA 

US = U t. ALPHA*D 


5 


DO 5 K=1 , NC 
DO 5 U=1,NPA 

US{ J,K)=U(U,K)+ALPHA*D( J,K) 


COMPUTE NEW TRAJECTORY AND PERFORMANCE INDEX 


CALL T.RAJTXS , US ) 

CALL PERFM (XS, US, PER) 


i'. ?N5B 
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IF(KOUT.GT.l) GO TO & 
WRI T£ ( 6 » 13 ) PER 
6 CONTINUE 


COMPUTE NEW GRADIENTS 
CALL ADJNT (XS i US t GS ) 
COMPUTE BETA 


7 


BN=0 » 

BD=0 c 

DO 7 K=lfNC 
DO 7 U=ltNPA 
BN=BN+GS ( J i K ) *GS ( J « K ) 
BD=BD+G( J,K)*G( J,K) 

beta=bn/bd 


ASSIGN NEW DIRECTION FOR LINE SEARCH 

DO 8-‘K=l»NC 
DO 8 J=1.NPA 

8 D(J,K)=-GS( J,K)+BET'A*D( JiK) 

MORE THAN MAXIMUM NUMBER OF TOTAL ITERATIONS 
IF(NOUT.LT, INDEX) GO TO 11 

MORE THAN MAXIMUM NUMBER OF INNER ITERATIONS 
IF(NOUT.LT.NPA) GO TO 9 . 


DOES PERFORMANCE INDEX IMPROVE ? 

IF ( PI EST-PER • LT. EP) GO TO 11 

PTEST-PER 

GO TO- -2 

MOVE NEW GRADIENT INT.O OLD GRADIENT 

9 DO 10 K=1,NC 
DO 10 J=1,NPA 
10 G ( J t K ) =GS ( K,) 

GO TO 4 


PRINT AND PLOT OPTIMAL TRAJECTORY 

11 KOUT=0 

CALL PRIN{XStUS.GSiPER«INDEX) 

DO 12 J=1 « NPA 

12 ALlSTf J)=FLOAT( J) 

CALL PLOTA(2.ALIST,US,NPAiNPA,NC) 
N=NPA+1 

„ CALL PL0TA(2,ALIST,XS.N,N,NX) 

13 FORMAT< ’+’i65X» ’PERFORMANCE INDEX IS 
STOP 

END 


*-i F15 . 4 ) 
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OIHENSION E XU0lIiflu(200,31 , G ( 200 ,3 ) , Y < 4 > , YN (4) , DLDX < 4 > , DFX ( 4 , 4 ) , 


? C0HM0N 3 PC0N(4) tTI4l , DELTA i KOUT , IND, NC , NX i NPA 
N=NPA+1 

INITIALIZE Y t K ) 

1 Y?K)=2.*PC0N IK) * (X(NtK ! -T IK) > AFLOAT CM) 

calculate the adjoint system in reverse time 

DO 2 J=ltNPA 
K = NPA+1- J 

define all derivatives of the form 

OFX ( NX i NX ) AND DFU(NX,NC) 

following are derivatives for the SECOND order jet engine model 

DEN=SQRT [ U ( K 1 1 ) **2+0 . 41688*X ( K i 2 ) **2-0 . 0899*U ( K • 1 ) *X ( K i 2 ) ) 

W=1.3009*X(K,2)-0.i3982S*U(K,l>-0.13982*DEl\j 

DWDU=-0. 139825-0. 06691* (2. *U(.K*1) -0 . 0899*X ( K t 2 ) >/DEN 
DWDX2=1. 3009-0. 06691* (0.63376*X<Ki2)-0.0899*U(K»l) i/DEN 
DFX {1 i 1 )=0 . 

DFU ( 2« 1 ) =ll258*?2?*U?Kl ll/ {X(Kil)*X(Kt2> ) -X ( K i 2 ) *DWDU ) 

CALCULATE DLX(NX) 

3 DLDX ( iT=2^*PCQN ( I ) * ( X ( K « I ) -T ( I ) ) *FLOAT < NPA )'* ( FLOA ! T ('K ) /FLOAT ( NPA ) ** 
79) 

CALCULATE PRESENT GRADIENT 

DO 4 I=liNC 
G ( K i I ) =0 • 

4 G?K?I ^gIki I)+Y(M)*DFU(M,I)*D£LTA 

CALCULATE NEXT Y(K) 

DO 5 M=1 * NX 
YN(M)=Y(H) 

DO 5 1=1 1 NX _ 

5 YN(M)=YN(M)+Y(I)*DFX(I»M)*DELTA 
DO 6 1=1 i NX 

YN( I )=YN( I )+DLDX(I)+Y(I ) 

6 YC I >=YN ( 1 » 

2 CONTINUE 

RETURN 

END 





Subroutine ADJNT computes the gradient arrays of a 
nominal trajectory <> The theory of the adjoint system 
equations is presented in section The user is 

required to calculate and define in &DJNT all derivatives 
of the form 

DFX ( J , K ) J=lj2, ...,NX K=1 # 2,...,NX 

DFU(J,K) Js=l t 2, »«*| NX K=1,2,...,NC 

where 

DPX( J,K) = — L • DPU( J,K) = L 

i** 5 u k 

The derivatives must be- expressed in the form which is 
illustrated by the jet engine problem example. The portion 
of ADJNT which ‘calculates the adjoint system and the gradient 
is completely general and need not be altered by the user. 

INPUTS : 

U - the control arrays 

X - the state variable arrays 

OUTPUT: 

G » the resulting gradient arrays 



SUBROUTINE CNSTR(X.UTRYiU) 

DIMENSION X { 4 ) »UTRY(3)»U(3) 

COMMON PC0N<4) , T (,H > , DELTA . KOUT« IND i NC i NX, NPA 
IFCUTRYCL) .GT.1.25*X<1) ) UTRYU)=X(1) 
IF(UTRYU) . GT » 1 « 25*X C 2 ) ) UTRY ( 1 ) =X ( 2 ) 
IFCUTRY(l) .LT.0.5) UTRY(1)=0.5 
U ( 1 ) =UTR Y ( 1 ) 

RETURN 

END 


Subroutine CNSTR tests a state variable n-tuple to 
determine the legality of the control, which is- reset at 
the nearest legal boundary if the control is- in an - 
unallowable region. 

This is a -specialized subroutine which implements the 
feasible directions modification. The user simply states 
his linear constraints in a similar manner. CNSTR is 
called only from the subroutine TRAJ. 


'INPUTS: 

UTRY - the present controls to be tes.te.d 
X - the present state variable n-tuple 
OUTPUT: 

U - the resulting allowable controls 



SUBROUTINE HOALF l X » U t D » ALP . H) „ „ ,, 

DIMENSION XC201., 4) , U ( 20 0 , 3 ) i D ( 20 0 1 3 ) ,UTRY{200,3> 
COMMON PC0N{4» tT(H) , DELI A tKOUT . INDi NC i NX i NPA 
DO 1 K=1 * NC 
DO 1 J=1 1 NPA 

UTRY(U.K)=U( J.K)+ALP*D< J,K) 

CALL TRAU(X,UTRY> 

CALL PERFM ( X i UTRY i H ) 

IF ( KOUT .01.1) BO TO 2 
IF ( IND.GT ,5) GO TO 2 
WRITE(b,3) ALP * H 
CONTINUE 

FORMAT (SX, 'ALPHA = 'iE15.5.' H( ALPHA) = ’iFIS.'f) 

RETURN 

END 


Subroutine HOALF calculates h(-a) = J(u + ad). It is 
called from LNSCH to determine the performance index 
resulting from possible control arrays. This subroutine 
is completely general and independent of the particular 
problem to be solved. It need' not be -altered by the user. 


INPUTS-: 

U - the nominal -cont"ol array 

X - the nominal state variable- arrays 

D - the direction along, which the line search Is made 

ALP - the present distance tested along D 

OUTPUT: 

H - the resulting performance index 
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WHICH UTILIZES A CUBIC POLYNOMIAL 

I T T°li E M®^ rc-iBBRKHiS OBTAINED FROM 
calculating the parameters needed FOR THE CUdlC FIT 

DIMENSION X (20114) ,U(200,3) , D ( 200 , 3) , G < 200 , 3 ) 

COMMON PCON (4) , T (4) , DELTA iKOUT , 1ND ,NC ,NX < NPA 
IND=IND+J. 

DETERMINE THE DERIVATIVE t H * ( 0 ) ) AND THE INITIAL • A-' GUESS 

HPO=0. 

AN=0. 

AD=0. 

DO 1 K=liNC 
DO 1 J-l i NPA 
HPO=HPO+G(U,K)*DC J,K) 

AN=AN+U( J,K)*U( J,K> 
ad=ao+0(U,K)*d: JiK) 

A=SQRT ( AN/'AD ) /10 . 

HO=PER 

CONTINUE 

CALL HOALF(XiUiDiAiHA) 

TEST THE RELATIONSHIP 8ETWEEN H(A) AND H(0) 

IF(HO.GT.HA) GO TO 4 

MAKE A BETTER GUESS FOR 'A* 


A=A/10 « 

IF(A.LT.1,E-12) GO TO 3 
GO TO 2 

3 ALPHA=0. 

GO TO a 

4 B=5.*A 

5 CONTINUE 

CALL HOALF(XiUiDiBiHB) 

TEST FOR RELATIONSHIP BETWEEN H C B ) AND HUM, _ALSO H ( B ) AND H ( A ) 


IF(HB.GT.HO) GO TO 6 
IF(HB+0.1.GT.HA) GO TO 7 


MOVE 'B* INTO ’ A ' AND MAKE A BETTER GUESS FOR *6* 


A=B 
HA=HB 
B=2.*B 
GO TO 5 


IF *B* IS A REASONABLE VALUE, USE THE CUBIC FIT TECHNIQUE 
TO FIND THE MINIMUM ALPHA 


AC=hA-HB*Ia/B )**3+H0* { ( A/B-) **3-1 •-) +HPO* ( A* ( A/B ) *#2- 
AC=AC/ { A**2-A*#3/B ) 

BC=(Ha-H0-HP0*A-AC*A*A)/A**3 

ALPHA= ( SbtRT ( rtC*AC-5 . *BC*HPO ) -AC ) /-( 3 . *BC ) 


•A) 


TEST TO DETERMINE IF PERFORMANCE INDEX IS TMFROVED BY ALPHA 


CALL HOALF(XiUiD,ALPHA,HS) 

if(hs.'lt.ha) go to a 

7 ALPhAsA 
6 CONTINUE 

PRINT RESULTS OF LNSCH DURING THE FIRST OUTER ITERATION 


IF(KOUT.GT.l) GO TO 9 
WRITE(failO) INO 
WRITE ( 6 1 11 ) ALPHA, HPO 
9 CONTINUE 

10 FORMAT (5X, • IND= ’,14) 

11 FORMAT (5X, 'THE CHOSEN 
RETURN 

END 


ALPHA IS ' , F10 • 9 , • HPO 


El5. 5 ) 
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SUBROUTINE PERFM ( X , U . PER ) 

DIMENSION X(201,4) ,U(200,3> 

COMMON PCON ( 4 ) ,T ( 4 ) , DELTA , KOUT i IND » NC t NX , NPA 
WGHT<K<N)=FLOAT(N) *{ FLOATt K > /FLOAT ( N) >**9 
N=NPA+1 
PER=0. 

DO 1 J=1 i N 
PEN=0. 

DO 2 K=1 i NX 

PEN=P£N+PCON (K)*(X(J»K)-TCK) ) **2 
PEN=PEN*WGHT( J»K) 

PER=PtR+PEN 

RETURN 

END 


Subroutine PERFM calculates the performance index for 
state variable and control arrays. This subroutine is 
completely general in nature and need not be altered by 
the user. 

INPUTS: 

U - the control arrays 

X - the state variable arrays 

OUTPUT: 

PER - the resulting performance index 


SSS33S 



SUBROUTINE PRI N ( X . U , G » PER . INDEX ) 

DIMENSION X(201.4),U(200.3).G(2G0.3) 

COMMON PCON14) .1(4) . DELTA i KOUT t IND i NC . NX . NP A 
WRITE ( 6 i 2 ) PER . INDEX 
IF(KOUT.GT.l) GO TO 4 
DO 1 U=1.NPA 

1 WRITE (6i 3)U. (X{ J.K) «K=1.NX) . (Ut J.IO .K=1.NX) , 
4 CONTINUE 

2 F0RMAK/.5X. 'PERFORMANCE INDEX IS »,F15.5»« 

3 FORMAT (HO . 10F10.4) 

RETURN 

END 


( G ( •« « K ) . K=1 . NX) 
ITERATION NUMBER '.15 


Subroutine PRIN prints the state variable and control 
arrays, and also any other desired- information. 

INPUTS: 

U - the control arrays 
X - the state variable arrays 
G - the gradient arrays 
PER - the performance index 

INDEX - the total number of iterations thus far 
OUTPUT: 

literal computer printout 


ORIGINAL PAGE IS 
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SUBROUTINE TRAJ(X.U-) 

DIMENSION X(201,4) ,U(-200,3) .UA(3) ,UT<3) »XT(4) 

COMMON PCON(H) »Tt4) .DELTA . KOUT i IND i NC i NX i NPA 

ENTER THE SYSTEM EQUATIONS HERE IN THE FORM 

Fl(XX»,..»XNXtUl tUNC) 

FNX ( xi i . . . » XNX . U1 1 . . . . UNC ) 

THE FOLLOWING EQUATIONS ARE EXAMPLES FROM THE 
SECOND ORDER JET ENGINE MODEL 

Fit XI ,X2.tUl»W)=37,78*W-38.448*Ul + 0.66849 
F2(XliX2*UltW>=1.258*(,Ul*Ul/(Xl*X2)-W*X2) 

VI ( U1 ,X2) =1.30 0 9*X2-l) . r39825*Ul 

V2(UliX2)=SG'RT<Ul*Ul+0.41686*X2*X2-0.069.9*X2*Ul) 

DO 1 J=1 1 NPA 
K=J+1 

CHOOSE STATE VARIABLE N-TUPLES AND CONTROLS. TO BE TESTED 
FOR THE FEASIBLE DIRECTIONS MODIFICATION 


2 

5 


4 


DO 2 M=1 * NX 
_XT(M)=X< J,M) 

DO 3 M=l‘iNC 
UT(M)=Ut J,M) 

CALL CNSTR ( XT « UT . UA ) 

THE RESULTING ALLOWABLE CONTROLS' 

DO 4 M=1»NC 
U ( J t M ) =UA (M) 

EULER INTEGRATION 


W=VltU( Jtl) tX{ J t '2) J-0V-l I 398'2*V2('U(UTI)-»-X4-JV2) ) 
x t k 1 1 ) =x (.j a > +F,X:(-X t j a% t ;x-(.J l *:2:)<vU.( J.«:i ),. iw-s-delta 
X t K 1 2 >'=X ( JIi-2 ) +F2 t'X C J vlhw’tft J,»,2.HTJ.( J »vi )!t Vn^DEE-W 
1 CONTINUE 
RETURN 
END 


Subroutine TRAJ computes- the state variable trajectory 
resulting from the control arrays; and’ a set of .initial 
conditions. The user enters the non-linear differential 
equations which describe the motion of the- system. 

Subroutine CNSTR is- called to* insure that the controls 
lie in an allowable region. 


INPUTS: 

U - the nominal control arrays 
X(l,) - a set of initial conditions 
OUTPUT: 

U - the allowable control arrays 
X - the resulting state variable trajectory 


OF' POOR QUALITY 
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